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ABSTRACT 

Many statistical methods have been proposed in the last years for analyzing the spatial 
distribution of galaxies. Very few of them, however, can handle properly the border 
effects of complex observational sample volumes. In this paper, we first show how to 
calculate the Minkowski Functionals (MF) taking into account these border effects. 
Then we present a multiscale extension of the MF which gives us more information 
about how the galaxies are spatially distributed. A range of examples using Gaussian 
random fields illustrate the results. Finally we have applied the Multiscale Minkowski 
Functionals (MMF) to the 2dF Galaxy Redshift Survey data. The MMF clearly in- 
dicates an evolution of morphology with scale. We also compare the 2dF real catalog 
with mock catalogs and found that ACDM simulations roughly fit the data, except at 
the finest scale. 

Key words: methods: data analysis - methods: statistical - large-scale structure of 
universe 



1 INTRODUCTION 

One of the main tenets of the present inflationary paradigm 
is the assumption of Gaussianity for the primordial density 
perturbations. This postulate forms the basis of present the- 
ories of formation and evolution of large-scale structure in 
the universe, and of its subsequent analysis. But it remains 
a hypothesis that needs to be checked. 

The most straightforward way to do that would be 
to f ollow t he de finition of Gaussian random fields (see, 
e.g, lAdled il98ll) l - their one-point probability distribu- 
tion and all many-point joint probability distributions of 
field amplitudes have to be Gaussian. This is clearly a too 
formidable task. Another way is to check the relationships 
between the correlation functions and power spectra of dif- 
ferent orders, which are well-defined for Gaussian random 
fi elds. This approach is fr equently used (see, e.g., a review 
in lMartfnez fc Saarl i2002l) '). A third method is to study the 
morphology of the cosmological (density) fields. One ap- 
proach to the morphological description relies on the so- 
called Minkowski functionals and is complementary to the 
moment-based methods because these functionals depend 
on moments of all orders. This procedure has been usually 
referred to as topological analysis. It has a qui te long his- 
tory a lready, starting with the seminal paper bv lGott et alJ 
( 1986), that deals with the genus, a quantity closely related 
to one of the four Minkowski functionals. The approach in 
the present paper lies within this latter framework; we de- 
scribe it in detail in Sec. [5] 



There are two different possibilities to develop a mor- 
phological analysis of galaxy catalogues based on Minkowski 
functionals. First, we can dress all points (galaxies) with 
spheres of a given radius, and study the morphology of 
the surface that is generated by the convex union of these 
spheres, as a function of the radius, which acts here as the 
diagnostic parameter. An appropriate theoretical model to 
compare with in this case is a Poisson point process. On 
the other hand, if we wish to study the morphology of the 
underlying realization of a random field, we have to restore 
the (density) field first, to choose an isodensity surface cor- 
responding to a given density threshold and to calculate 
its morphological descriptors. In this approach the density 
threshold (or a related quantity) acts as the diagnostic pa- 
rameter. The theoretical reference model is that of a Gaus- 
sian random field, and the crucial point here is to properly 
choose a restoration method that provides a smoothed un- 
derlying density field that should be fairly sampled by the 
observed discrete point distribution. 

Starting from the original paper on topology by 
iGott et all (1986), this task has been typically done by 
smoothing the point distribution with a Gaussian kernel. 
The choice of the optimal width of this kernel has been 
widely discussed; it is usually taken close to the correlation 
length of the point distribution. However, Gaussian smooth- 
ing is not the best c hoice for morp hological studies. As we 
have shown recently iMartfnez et al.l2005l) . it tends to intro- 
duce additional Gaussian features even for manifestly non- 
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Gaussian density distributions. Minkowski functionals are 
very sensitive to small density variations, and the wings of 
Gaussian kernels could be wide enough to generate a small- 
amplitude Gaussian ripple that is added to the true density 
distribution. Such effect could be alleviated by using com- 
pact adaptive smoothing kernels, as we show below. 

It is well known that large-scale cosmological fields have 
a multi-scale structure. A good example is the density field; 
it includes components that vary on widely different sales. 
The amplitudes of these components can be characterized by 
the power spectrum; the present determinations encompass 
the frequency interval 0.01-0.8 /i/Mpc, 1 which corresponds 
to the scale range from 8 to over 600 Mpc/h. Fig. □ shows 
the power spectrum for the 2dF galaxy redshift survey (2dF- 
GRS). 

As cosmological densities have many scales and widely 
varying amplitudes, density restoration should be adaptive. 
Different methods exist to adaptively smooth point distri- 
butions to estimate from them th e underlying density field. 
ISchaap fc van de WeveaerU J2(3Q0) have introduced the De- 
launay Tessellation Field Estimator (DTFS) which adapts 
itself to the point configuration even when anisotropies 
are present. The method starts by considering the Delau- 
nay tessellation of the point process, then we can estimate 
the density at those points using the contiguous Voronoi 
cells, and finally, we should interpolate to obtain the den- 
sity in the whole volume. Intricate point patterns have 
been successfully smoothed using this method and applica- 
tions to particle hydrodynamics provide good per formance 
iPelupessv et alJl2003l) . lAscasibar fc Binnevl i2005l) have re- 
cently introduced a novel technique based on a different 
partition of the embedding space. These authors use mul- 
tidimensional binary trees to make the partition and latter 
apply adaptive kernels within the resulting cells. Finally, it 
is well known that wavelets provide a l ocalized (compact- 
kerne l) adaptive restoration method iStarck fc Murtagh 
2 QQ2 ). We have applied, in a previous paper jM^irtrnez^t^l 



2005), a wavelet based denoising technique to the 2dFGRS. 



As a result we found that the morphology of the galaxy 
density distribution in the survey volume does not follow a 
Gaussian pattern, in contrast to the usual results in which 
deviations of Gauss ianity are not clearly detected (see, e.g. , 
iHovle et al.1 (l2002ft for the 2dFGRS and lPark etaH fcOOBft 
for the Sloan Digital Sky Survey (SDSS)). 

By the way, adaptive density restoration methods are 
probably the best for calculating partial Minkowski func- 
tionals, to describe the morphology of singl e large-scale den- 
sity en hancements (superclusters; see, e.g.. IShandarin et alJ 
(2004)). Partial functionals can be used to characterize the 
inner structure (d umpiness) and shapes (via shapefinders 
iSahni et aljH99 8)) of superclusters. As Minkowski func- 
tionals are additive, partial functionals can be, in principle, 
combined to obtain global Miknkowski functionals for the 
whole catalogue volume. But if we want to check for non- 
Gaussianity, direct calculation of global Minkowski function- 
als is more simple and straightforward. When combining 
partial functionals, estimating the mean densities and vol- 
ume distributions for the full sample is a difficult problem. 



1 As usual, h is the present Hubble parameter, measured in units 
of 100 km/sec/Mpc. 
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Figure 1. The matter density power spectrum for the 2dF GRS. 
Data courtesy of W. Percival and the 2dF GRS team. 



Now, although a single adaptively found density dis- 
tribution represents the cosmological density field better, it 
could not be the best tool for comparing theories with obser- 
vations. Theories of evolution of structure predict that Gaus- 
sianity of the original density distribution is distorted during 
evolution, and this distortion is scale-dependent. Accord- 
ing to the present paradigm, evolution of structure should 
proceed with different pace at different scales. At smaller 
scales signatures of gravitational dynamics should be seen, 
and traces of initial conditions could be discovered at larger 
scales. In a single density field, containing contributions from 
all scales, these effects are mixed. Thus, a natural way to 
study cosmological density fields is the multi-scale approach, 
scale by scale. This has been don e in the past by u sing a se- 
ries of kernels of different width llPark e t al. 2005), but this 
method retains considerable low-frequency overlap. A better 
way is to decompose the density field into different frequency 
(scale) subbands, and to study each subband separately. 

The simplest idea of separation of scales by using differ- 
ent Fourier modes does not work well, at least for morpho- 
logical studies (studies of shapes and texture). Describing 
texture requires knowledge of positions, but Fourier modes 
do not have positions, their position is the whole sample 
space. A similar weakness, to a smaller extent, is shared by 
discrete orthogonal wavelet expansions - their localisation 
properties are better, but vary with scale, and large-scale 
modes remain badly localised. 

This leads to the conclusion that the natural candi- 
dates for scale separation are shift-invariant wavelet systems, 
where wavelet amplitudes of all scales are calculated for each 
point of the coordinate grid. These wavelet decompositions 
are redundant - each subband has the same data volume as 
the original data. 

For such a scheme, direct calculation of low- frequency 
subbands would require convolution with wide wavelet pro- 
files, that could be numerically expensive. A way out is the 
d trous (with holes) trick, where convolution kernels for the 
next dyadic scale are obtained by inserting zeros between 
the elements of the original kernel. In this way, the total 
number of non-zero elements is the same for all kernels, and 
all wavelet transforms are equally fast. 

Now, as usual with wavelets, there is considerable free- 
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dom in choosing the wavelet kernel. A particularly useful 
choise is the wavelet based on a B$ spline scaling function 
(see Appendix A) . In this case, the original data can be re- 
constructed as a simnple sum of subbands, without extra 
weights, so the subband decomposition is the most natural. 

In the present paper we combine the a trous representa- 
tion of density fields with a grid-based algorithm to calculate 
the Minkowski functionals (MF for short), and apply it to 
the 2dFGRS data. 

Section [5] describes how to compute the MF for com- 
plex data volumes and how to extend such an approach in a 
multiscale framework using the wavelet transform. Sectional 
describes our observational data while section 2] evaluates 
the multiscale MF on Gaussian random field realizations for 
which the analytical results are known. Section |5] presents 
our results for the 2dFGRS data. These results are com- 
pared with the multiscale MF calculated for 22 mock sur- 
veys and about 100 Monte-Carlo simulations of Gaussian 
random fields. We list the conclusions in section [5] 



2 MORPHOLOGICAL ANALYSIS 
2.1 Definition 

An elegant description of morphological characteristics 
of density fields is given by Minkowski functionals 
iMecke et al.lll993) . These functionals provide a complete 
family of morphological measures. In fact all additive, mo- 
tion invariant and conditionally continuous 2 functionals de- 
fined for any hypersurface are linear combinations of its 
Minkowski functionals. 

The Minkowski functionals describe the morp hology of 
iso-density surfaces jMinkowskll9od : lTomita ^?90|), and de- 
pend thus on the specific density level (see ISheth fc Sahnil 
(2005) for a recent review). Of course, when the original 
data are galaxy positions, the procedure chosen to calcu- 
late densities (smoot hing) will also determine the result 
iMartfnez et al.ll2005l) . Generally, convolution of the data 
with a Gaussian kernel is applied to obtain a continuous 
density field from the point distribution. An alternative ap- 
proach starts from the point field, decorating the points with 
spheres of the same r adius, and studying the morphology of 
the re sulting surface dSchmalzing et al . 1996: lKerscher et alJ 
Il997f) . This approach does not refer to a density; we cannot 
use that for the present study. 

The Minkowski functionals are defined as follows. Con- 
sider an excursion set F<j> of a field 0(x) in 3-D (the set of 
all points where 0(x ^ <f>o). Then, the first Minkowski func- 
tional (the volume functional) is the volume of the excursion 
set: 



Vi(4>o) 



VoW>o) 



The second MF is proportional to the surface area of the 
boundary 5Fa of the excursion set: 



2 The functionals are required to be continuous only for compact 
convex sets; we can always represent any hypersurface as unions 
of such sets. 
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The third MF is proportional to the integrated mean curva- 
ture of the boundary: 
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where Ri and R2 are the principal curvatures of the bound- 
ary. The fourth Minkowski functional is proportional to the 
integrated Gaussian curvature (the Euler characteristic) of 
the boundary: 
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The last MF is simply related to other known morphological 
quantities 

Vs = X=i(l-G), 

where \ is the Euler characteristic and G is the topological 
genus, widely used in the past study of cosmological density 
distributions. The functional V3 is a bit more comfortable 
to use - it is additive, while G is not, and it gives just twice 
the number of isolated balls (or holes). Instead of the func- 
tionals, their spatial densities Vi are frequently used: 

n{f) = Vi{f)/v, i = o,...,3, 

where V is the total sample volume. The densities allow us 
to compare the morphology of different data samples. 

The original argument of the functionals, the density 
level po, can have different amplitudes for different fields, 
and the functionals are difficult to compare. Because of that, 
normalised arguments are usually used; the simplest one is 
the volume fraction f v , the ratio of the volume of the ex- 
cursion set to the total volume of the region where the den- 
sity is defined. Another, similar argument is the mass ratio 
fm, which is very useful for real, positive density fields, but 
is cumbersome to apply for realizations of Gaussian fields, 
where the density may be negative. The most widely used 
argument is the Gaussianized volume fraction v, defined as 
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exp(-t /2) dt. 
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For a Gaussian random field, v is the density deviation from 
the mean, divided by the st andard deviation. T his argument 
was introduced already bv lOott et alJ (jlflSfiT )') . in order to 
eliminate the first trivial effect of gravitational clustering, 
the deviation of the 1-point pdf from the (supposedly) Gaus- 
sian initial pdf. Notice that using this argument, the first 
Minkowski functional is trivially Gaussian by definition. 

All the Minkowski functionals have analytic expres- 
sions for iso-density slices of realizations of G aussian random 
fields. For three-dimensional space they are <|Tomital H990): 
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where <&(•) is the Gaussian error integral, and A is deter- 
mined by the correlation function £(r) of the field: 
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2.2 Numerical algorithms 



Several algorithms are used to calculate the Minkowski func- 
tionals for a given density field and a given density threshold. 
We can either try to follow exactly the g eometry of the iso - 
density surface, e.g., using triangulation feheth et alJl2003» . 
or to approximate the excursion set on a simpl e cubic lat- 
tice. T he algorithm that was proposed first by iGott et alJ 
(1986), uses a decomposition of the field int o filled and 
empt y cells, and another popular algorithm dColes et alJ 
1996) uses a grid-valued density distribution. The lattice- 
based algorithms are simpler and faster, but not as accurate 
as the triangulation codes. 

We use a simple grid-based algorithm, that makes use 
of integral geometry ( Crofto n's intersection formula, see 
ISchmalzing fc Buchertl Jl997)). We find the density thresh- 
olds for given filling fractions by sorting the grid densities, 
first. Vertices with higher densities than the threshold form 
the excursion set. This set is characterised by its basic sets of 
different dimensions - points (vertices) , edges formed by two 
neighbouring points, squares (faces) formed by four edges, 
and cubes formed by six faces. The algorithm counts the 
numbers of elements of all basic sets, and finds the values of 
the Minkowski functionals as 
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where a is the grid step, / is the filling factor, iVo is the 
number of vertices, Ni is the number of edges, N2 is the 
number of squares (faces), and N3 is the number of basic 
cubes in the excursion set for a given filling factor (density 
threshold) . The formu l a (171 was first used in cosmological 
studies bv lColes etafl l|l996h . 



2.3 Biases 

The algorithm described above is simple to program, and is 
very fast, allowing the use of Monte-Carlo simulations for 
error estimation. 

However, it suffers from discreteness errors, which are 
not large, but annoying, nevertheless. An example of that is 
given in Fig. |5J where we show the Vi functional, calculated 
by the above recipes for a periodic realization of a Gaussian 
field (the dashed line). As we see, it has a constant shift in v 
over the whole range. This shift is due to the fact that when 
we approximate iso-density surfaces by a discrete grid, the 
vertices that compose the surface lie in a range of densities 
starting from the nominal one. This effect can easily be cal- 
culated because this bias will show up as a constant shift in 
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Figure 2. The V\ functional for a realization of a Gaussian ran- 
dom density field in a periodic 256 3 cube (upper panel). Full line 
shows the theoretical prediction for this realization, dotted line - 
the standard one-excursion-set estimate, and dashed line - the av- 
erage of the functional over two excursion sets, The lower panel 
shows the difference between the estimates and the theoretical 
prediction; the lines encode the estimates as above. 



v for a Gaussian density field, as observed. Other function- 
als (V2 and V3) suffer similar shifts, with smaller amplitudes, 
and these are not easy to explain. 

There is, fortunately, another and simple possibility to 
fight these errors. The standard way is to approximate an 
iso-density surface by the collection of vertices that have 
densities p ^ pi, where pi is the threshold density. But an- 
other surface, formed by the vertices with p < pi, is as good 
an approximation to the iso-density surface as the first one. 
Thus, the natural way to calculate the Minkowski function- 
als is to run the algorithm twice, swapping the marks for 
the excursion set, and averaging the values of the function- 
als obtained. The last step is justified, as the Minkowski 
functionals are additive. The averaging rules are: 



V 

v 2 



(y«+y ot -vf>) 1 
(Vw + vf>) A 
(y (1) -y (2) ) 1% 
(y™ + y 3 (2) ) A 



Here the upper indices (1) and (2) denote the original and 
complementary excursion sets, respectively, and Vt ot is the 
total number of grid cubes in the data brick. The minus sign 
in the formula for the third functional (V2) accounts for the 
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fact that the curvature of the second surface is opposite to 
that of the first one. 

The Vi functional calculated this way is shown in Fig.|3 
by the full line; we can compare it with the theoretical pre- 
diction for Gaussian fields (©, the dotted line). The coinci- 
dence of the two curves is very good; the only slight deviation 
is at v w 0, where the Gaussian surface is more complex. We 
have to stress that the Gaussian curve is not a fit; the pa- 
rameter A that determines the amplitude of the curve was 
found directly from the data, using the relations £(0) = {p 2 } 
and £"(0) = {p 2 i), where £(r) is the correlation function and 
p t i is the derivative of density at a grid vertex in one of 
the coordinate directions. The good match of these curves 
shows also that the Gaussian realization is good, which is 
not simple to model. The averaging works as well for the two 
other functionals; this is shown in Fig. [3] There are slight 
deviations from the theoretical curve for V2 around v « 1 
and for V3 at v w 0; these may be intrinsic to the particular 
realization, as the number of 'resolution details' diminishes 
when the order of the functional increases. 

We also see that the higher the order, the closer are the 
one- and two- excursion set estimates. So, even if we are in- 
terested only in the topology of the density iso-surfaces, we 
should correct for the border effects, all Minkowski function- 
als are used, and it is important that they were unbiased. 

There is a natural restriction on the grid steps - the 
grid has to be fine enough to resolve the details of the 
density field. The previous figures (Figs. EH show the 
Minkowski functionals obtained for the case of the Gaus- 
sian field smoothed by a Gaussian filter with a = 3 grid 
steps, being the smoothing radius R ~ 2a = 6. 

2.4 Border corrections 

As we have seen above, we can obtain good estimates of 
the Minkowski functionals for periodic fields. The real data, 
however, is always spatially limited, and the limiting sur- 
faces cut the iso-density surface. An extremely valuable 
property of Minkowski functionals is that such cuts can be 
corrected for. Let us assume that the data region (window or 
mask) is big enough relative to the typical size of details, so 
that one can consider the fiel d inside the mask homoge neous 
and isotropic. For this case, ISchmalzing et all il996l) show 
that the observed Minkowski functionals for the masked iso- 
surface Mi(DCiW) can be expressed as a combination of the 
true functionals Mi(D) and those of the mask Mi(W): 

Mi{D n W) = i ( ) ) M > ( D ) M ^ W> (8) 

3=0 ^ ' ' 

where V is the total volume inside the mask. Note that the 
functionals Mj d i ffer fr om the usual Vi by normalisation. 
ISchmalzing et all (I1996T) derive the relation iJHJ for a collec- 
tion of balls. Here we have applied it to iso-density surfaces 
and for the true values of the functionals, we get 

Mj(D) _ MdpnW) jp/ i \ Mj{D) Mj-jjW) 
V M (W) j^ \j J v M (W) ' ( ' 

The relation between Mj-s and the usual Vi-s is 
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Figure 3. The V2 (upper panels) and V3 (lower panels) func- 
tionals for a realization of a Gaussian random density field in a 
periodic 256 3 cube. In the larger panels full lines show the theo- 
retical predictions for this realization, dotted lines - the standard 
one-excursion-set estimates, and dashed lines - the averages of the 
functionals over two excursion sets, The smaller panels show the 
differences between the estimates and the theoretical predictions; 
the lines encode the estimates as above. 
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where uij is the volume of a j-dimensional unit ball, and d 
is the dimension of the space. For Minkowski functionals in 
three-dimensional space, the explicit relations are: 

M = V ; Afi = fvi; M 2 = ^-V 2 ; M 3 = ^-V 3 . 

4 Air Atv 

Using © and replacing Mj-s by Vi, we arrive at the following 
correction chain: 
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Here Vi(i/) denote the observed (raw) values of Minkowski 
functionals, and Vi{v) denote the corrected densities. 

We tested these corrections with our original Gaussian 
realization, masked at all faces. The correction for the sec- 
ond Minkowski functional wi is practically perfect. The cor- 
rected version of v 2 is also close to the original for all the 
argument range, and only a little higher than the original. 
The higher the order of the functional, the more difficult it 
is to correct for the borders, as small errors from the lower 
orders accumulate. The discrepancy with the corrected ver- 
sion of vz and the original estimates is the largest amongst 
the three densities, but it balances well the amplitudes of 
the maxima, and is only a little lower than the original for 
low densities [y ~ 1.5). There are practically no differences 
from the original at the high-density end. 

A note on the use of masks in practice: we ensure that 
there is at least one- vertex thick mask layer around our data 
brick. This allows us to assume periodic borders for the brick 
itself. And there is also another way to use the mask, ignor- 
ing the vertices in the mask, not building any elements from 
the vertices in the data region to the mask vertices. Then we 
do not have to apply the correction chain 11U> and do not 
have to build the basic sets in the mask region. The latter 
fact makes the algorithm about twice as fast for the 2dF data 
(the data region occupies only a fraction of the encompassing 
brick). We compared this version with the border-corrected 
algorithm described above (see Appendix B) and found that 
it gives slightly worse results for Gaussian realizations, so we 
dropped it. The present algorithm is fast enough, taking 12 
seconds for a iso-level for a 256 3 grid on a laptop with the 
Intel Celeron 1500 MHz processor. 

As the data masks are complex (see Fig. [SJ, we should 
test the border corrections for real masks, too. Fig. 01 shows 
the effect of border corrections as used with the data mask 
for the 2dF NGC sample volume (see section [3J. As the 
corrections give densities of functionals, we will show the 
densities from this point on. The density field in this volume 
was generated by simulating a Gaussian random field for a 
256 x 256 x 64 periodic brick, combining these bricks to cover 
the sample extent, and masking this realization with the 
Northern data mask. Smaller bricks were combined because 



the spatial extent of the data was too large for the available 
core memory to generate a single FFT brick to cover it; 
as the brick is periodic, the realization remains Gaussian. 
We show the raw densities of the Minkowski functionals as 
dashed lines, the corrected versions with full lines, and the 
densities for the original brick with dotted lines. 

We see that the density v\ of the V\ functional is re- 
stored well, apart from a slight deviation near v = 0. The 
density v 2 is also corrected well, only its maximum ampli- 
tude is slightly smaller than that for the original brick. The 
restoration is almost perfect for V3 , with the same small am- 
plitude problem than for the two other densities. The fact 
that the restoration works so well is really surprising - first, 
our mask is extremely complex, and, secondly, our realiza- 
tion of the Gaussian random field is certainly not exactly 
homogeneous and isotropic inside the sample volume. We 
generate realizations of random fields in this work by the 
FFT technique; these fields are homogeneous, but isotropic 
only for small scales, not for scales comparable to the brick 
size. 

The importance of the good restoration of the 
Minkowski functional means that when checking theoretical 
predictions, we can directly compare observational results 
with the predictions, and do not have to use costly Monte- 
Carlo simulations. 



2.5 Multiscale Minkowski Functionals 

A natural way to study cosmological fields is the multi-scale 
approach, scale by scale. According to the present paradigm, 
evolution of structure should proceed with different pace at 
different scales. At smaller scales signatures of gravitational 
dynamics should be seen, and traces of initial conditions 
could be discovered at larger scales. 

The matter density in the universe is formed by per- 
turbations at all scales. In the beginning they all grow at 
a similar rate, but soon this rate becomes scale-dependent, 
the smaller the scale, the faster it will go nonlinear and non- 
Gaussian. Thus it is interesting to decompose the density 
(and gravitational potential, and velocity) field into differ- 
ent scales and check their Gaussianity (and other interesting 
characteristics). 

Using the wavelet transform as it is described in 
Appendix A, we obtain a set of wavelet scales W = 
Wo,..., wj, cj+i, and each scale Wj(x, y, z) corresponds to the 
convolution product of the observed galaxies with a wavelet 
function ipj, where ipj(x,y, z) — ^(^-, 53-, and tp is the 
analyzing wavelet function described in Appendix A. Now, 
we can apply the MF calculation at each scale independently, 
and we get four MF values per scale using Eq. |7| The set 
(Vj,o, Vj,i> Vj',3) will denote the MF at scale j. Note that 
in this framework, we do not have to convolve the data any- 
more with a Gaussian kernel, avoiding the delicate choice of 
the size of the kernel bandwidth. 



3 THE DATA 

There exist two large-volume galaxy redshift surveys at the 
moment, t he Two Degree Fie ld Galaxy Redshift Survey 
(2dFGRS) JColless et alJl2003ft and the Sloan Digital Sky 



7 



6 | , , , , , , , 1 0.8 | , , , , , , , 1 1.5 




Figure 4. Demonstration of border corrections for complex borders. The raw densities of Minkowski functionals for the 2dFGRS NGP 
sample volume 252 3 , cut from a periodic 256 2 X 64 realization of a Gaussian random field (dotted lines) are shown together with the 
border-corrected estimates (dashed lines) and the estimates for the original brick (full lines). Upper panels show the densities, and smaller 
lower panels show the differences between the densities for the sample volume and those for the brick (for both the raw and corrected 
cases, and the same line types are used as in the upper panels). The densities of the second MF v\ are shown in the left panels, the 
densities of the third MF V2 — in the middle panels, and the densities of the fourth MF — in the right panels. 



Survey (SDSS) <York et alJ hoOOl) . The 2HFGRS is com- 
pleted; although the SDSS is not, its data volume has al- 
ready surpassed that of the 2dFGRS. We shall use in our 
paper the 2dFGRS dataset; it is easier to handle, and we 
make use of the mock catalogues created to estimate the 
cosmic variance of the data. 

The galaxies of the 2dFGRS h ave been selected fro m an 
earlier photometric APM survey (iMaddox et al]ll996t) and 
its extensions. The survey covers about 2000 deg 2 in the 
sky and consists of two separate regions, one in the North 
Galactic Cap (NGC) and the other in the South Galactic 
Cap (SGC), plus a number of small randomly located fields; 
we do not use the latter. The total number of galaxies in the 
survey is about 250,000. The depth of the survey is deter- 
mined by its limiting apparent magnitude, which was chosen 
to be bj = 19.45. Caused by varying observing conditions, 
however, this limit depends on the sky coordinates, varying 
almost the full magnitude. Another cause of non-uniformity 
of the catalogue is its spectroscopic incompleteness - as the 
fibres used to direct the light from a galaxy image in the 
focal plane to the spectrograph have a finite size, a num- 
ber of galaxies in close pairs were not observed for redshifts. 
However, these corrections can be estimated; the 2dFGRS 
team has made public the programs that calculate the com- 
pleteness factors and magnitude limit, given a line-of-sight 
direction. 

The 2dFGRS survey, as all redshift surveys, is 
magnitude-limited. This means that the density of observed 
galaxies decreases with distance; at large distances only in- 
trinsically brighter galaxies can be seen. For certain sta- 
tistical studies (luminosity functions, correlation functions, 
power spectra) this decrease can be corrected for. For tex- 
ture studies there are yet no appropriate correction meth- 
ods, and maybe, these do not exist, as the scales of the 
details that can be resolved are inevitably different, and 
small-scale information is certainly lost at large distances. 
The usual approach is to use volume-limited subsamples 



extracted from the survey. In order to create such a sam- 
ple, one chooses absolute magnitude limits, and retains only 
the galaxies with absolute magnitudes between these limits. 
This discards most of the data, but assures that the spatial 
resolution is the same throughout the survey volume (tak- 
ing also account of the possible luminosity evolution and 
K-correction). 

The 2dFGRS team has created such catalogues and 
used them to study higher-order corre lation functions 
jNorberg et alJ l2002t ICroton et ail l2004allbh . They kindly 
made these catalogues available to us, and these constitute 
our main data. These catalogues span one magnitude each, 
from M = -17 + 5 log 10(h) until M = -21 + 5 log 10(h). 
The catalogues for least bright galaxies span small spatial 
volumes, and those for the brightest galaxies are sparsely 
populated. Thus we chose for our work the catalogue span- 
ning the magnitude range —20 M — 5 log 10(h) —19, this 
is the most info r mative . This has been also the conclusion 
of lCroton et alJ J2004bl) . We shall call this sample 2dF19. 
This sample, as all 2dF volume-limited samples, consists of 
two spatially distinct subsamples, one in the North Galac- 
tic Cap region (2dfl9N) and another in the South Galactic 
Cap region (2dfl9S). The sample lies between 61.1 Mpc/h 
and 375.6 Mpc//i; the general features of the subsamples are 
listed in Table [T] 

In a previous pap er on the morphology of the 2dFGRS 
iMartmez et aljl2005l) we extracted bricks from the data to 
avoid the influence of border effects. This forced us to use 
only a fraction of the volume-limited samples. This time we 
tried to use all the available data, and succeeded with that. 

The 2dFGRS catalogue is composed of measurements in 
a large number of circular patches in the sky, and its foot- 
print in the sky is relatively complex (see the survey web- 
page http://www.mso.anu.au/2dFGRS). Furthermore, due 
to the variations in the final magnitude limit used, the cat- 
alogue depth is also a function of the direction. In order to 
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Table 1. The 2dF volume-limited catalogues used. 



sample 


galaxies 


ra limits (deg) 


dec limits 


Vol (10 6 


dmean 






(deg) 


(deg) 


Mpc 3 h- 3 ) 


(Mpc/ft) 


2dF19N 


19080 


147.0 223.0 


-6.4 2.6 


2.75 


5.24 


2dF19S 


25633 


-35.5 55.2 


-37.6 -22.4 


4.43 


5.57 




Figure 5. Survey masks for the 2dFGRS volume-limited sample 
2dF19. Upper panel — the Northern subsample, lower panel - 
the Southern subsample. Spatial orientations are chosen to better 
visualise the volumes. 

use all the data for wavelet and morphological analysis, we 
had to create a spatial mask, separating the sample volume 
from regions outside. We did it by creating first a spatial 
grid for a brick that surrounded the observed catalogue vol- 
ume, calculated the sky coordinates for all vertices of the 
brick, and used then the software provided by the 2dFGRS 
team to find the correction factors for these directions. The 
completeness factor told us if the direction was inside the 
sample footpath. If it was, we found the apparent magnitude 
of the brightest galaxies (with M = M m i n ) of our sample for 
that distance and checked, if it was lower than the sample 
limit. If it was, the grid point was included in the mask. For 
the last comparison we had to change from the comoving 
grid distance to the luminosity distance. We assumed the 
fiiW = 0.3, £Ia — 0.7 cosmology for that and interpolated a 
ta bulated relation betwe en the two distances. As explained 
in iNorberg et alJ <l2002tl , the 2dF volume- limited samples 
were built using the k + e-correction as dependent on the 
spectral type of a galaxy. We do not have such a quantity 
for the mask, so we tuned a little the bright absolute mag- 
nitude limit, checking that the mask should extend as far as 
the galaxy sample. For that, we had to increase the effective 
bright absolute magnitude limit by AM = 0.25. The nearby 
regions of the mask were cut off at the nearest distance limit 
for the observational sample. 



The original survey mask in the sky includes holes 
around bright stars; these holes generate narrow tunnels 
through our spatial mask. As the masks have a complex 
geometry, we did not want to add discreteness effects due 
to resolving such tunnels. We filled them in by counting the 
number of neighbours in a 3 3 cube around non-mask points 
and, if the neighbour number was larger than a chosen limit, 
assigning the points to the mask. We chose the required 
number n of neighbours to avoid filling in at flat mask bor- 
ders (n = 9 is enough for that) and iterated the procedure 
until the tunnels disappeared. This was de termined by vi- 
sual checks (using the 'ds9' fits file viewer lljove fc Mandell 
120031) ). 

We show the 3-D views of our masks in Fig. |5] The 
mask volumes are, in general, relatively thin curved slices 
with heavily corrugated outer walls. These corrugations are 
caused by the unobserved survey fields. Also, the outer edges 
of the mask are uneven, due to the variations in the survey 
magnitude limit. One 3-D view does not give a good im- 
pression of the mask; the slices we shall show below will 
complement these. 

As mentioned in Appendix A, in order to apply the 
wavelet convolution cascade, the initial density on the grid 
should be extirpolated by the chosen scaling function. We 
chose our initial grid step as 1 Mpc/h, and used the ker- 
nel for extirpolation. In order to have a better scale coverage, 
we repeated the analysis, using the grid step v2 Mpc/ft. The 
smoothing scale (the spatial extent of the kernel) is 4 grid 
units, the smoothing radius corresponds to 2 units. When 
calculating the densities, we used the spectroscopic com- 
pleteness corrections ej Bp , included in the 2dFGRS volume- 
limited catalogues, and weighted the galaxies by the factor 
w = 1/ c: . Most of the weights are close to unity, but a few 
of them are large. In order not to 'overweight' these galax- 
ies, we fixed the maximu m weight level as w = 2. The same 
procedure was chosen bv lCroton et alJ l|2004al) . 



4 GAUSSIAN FIELDS 

The filters used to perform wavelet expansion are linear, 
and thus should keep the morphological structure of Gaus- 
sian fields; the Minkowski functionals should be Gaussian 
for any wavelet order. This is certainly true for periodic 
densities, but for densities restricted to finite volumes the 
boundary conditions can introduce correlations. The most 
popular boundary condition - reflection at the boundary - 
will keep the density field mostly Gaussian for brick masks. 
Our adopted zero boundary condition will certainly work 
destroying Gaussianity, as the random field which is zero 
outside a given volume and has finite values inside is cer- 
tainly not Gaussian. 

We compare the effect of the mirror and zero boundary 
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Figure 6. Densities of the Minkowski functionals (from the second at the left until the fourth at the right) for the same Gaussian 
realization for the wavelet orders 3 and 4 (lower wavelet orders give higher amplitudes). The case for the wavelets generated using the 
mirror boundary conditions is shown by full lines, for zero boundary conditions dotted lines are used. The amplitudes of the functionals 
for the wavelet order 4 are rescaled by 2 for V2 and by 4 for to show more details. 
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Figure 7. Multi-scale decomposition of a realization of a Gaussian random field in the 2dfl9N volume mask, for the z = 34 Mpc//i slice 
(the data and the first orders). The upper panel column shows the scaling orders, with the original density at the right. The lower panel 
shows the wavelet orders, with the lowest order at the right. 



conditions in Fig. El for high wavelet orders that are yet not 
dominated by noise, calculated for a Gaussian realization in 
a 256 3 cube with the a = 3 Gaussian smoothing described 
above, and masked at all borders by a layer two vertices 
thick (thus the effective volume of the data cube is 252 3 ) 

Fig. shows that the second Minkowski functional is 
certainly better restored by the wavelets obtained by using 



mirror boundary conditions. These conditions lead to the 
functionals that are symmetric about v = 0, as it should be, 
while the functionals obtained by applying zero boundary 
conditions display a shift towards smaller values of v. How- 
ever, although our brick data should prefer mirror boundary 
conditions, it is not easy to say which boundary conditions 
give better estimates for the remaining two functionals. Both 
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Figure 8. Multi-scale decomposition of of a Gaussian random field in the 2dfl9N volume mask (continued). Upper panel - scaling 
orders, lower panel — wavelet orders, the highest orders at the left. The last scaling solution shows already strong effects of the boundary 
conditions. 



cases have comparable errors (look at the amplitudes at ex- 
trema, which should be equal). Hence, other considerations 
have to be used. As our data mask is complex, with corru- 
gated planes and sharp corners, the mirror boundary condi- 
tions will amplify the corner densities and propagate them 
inside, while the zero boundary conditions will gradually re- 
move the influence of the boundary data. Thus we have se- 
lected zero boundary conditions for the present study; these 
are probably natural for all observational samples. 

We have seen that the wavelet components of the multi- 
scale decomposition of a realization of a Gaussian random 
field remain practically Gaussian for simple sample bound- 
aries. This means that testing for Gaussianity is straight- 
forward. However, we have to assess the boundary effect for 
our application, where the boundaries are extremely com- 
plex. We shall demonstrate it on the example of the North- 
ern mask. For that, we generated a realization of a Gaus- 
sian random field for a volume encompassing the mask, 
as described in section and masked out the region out- 
side the NGP data volume. As we want to see the ef- 
fects that could show up in the data, we used the stan- 
dard dark matter power sp ectrum for the ACDM cosmology 
jKlvpin fc Holtzmarll997h . for the cosmological parameters 
fl = 0.3, n A = 0.7,n ba r = 0.026, h = 0.7 (this is pretty 



0.5 



-0.5 



-1.5 















"H... 8W2 


.{/ 




\ 2W1 












\W0 ./ 

l i *-T' 









-1.5 



-0.5 



0.5 



1.5 



Figure 9. The density 113 of the fourth MF for the wavelet decom- 
position of a realization of a Gaussian random field in the 2dfN19 
mask. The legends in the figure show the wavelet order and the 
scaling factor (2W1 denotes the functional for the wavelet order 1, 
multiplied by 2). The legend NG denotes the original realization. 
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Figure 10. Multi-scale decomposition of the 2dfl9N volume-limited sample, for the z = 39 Mpc/h slice, and for the grid step \f2 Mpc/h. 
The original density is shown at top right, the last scaling order at bottom left, the wavelet orders in between. The wavelet orders increase 
from right to left and from top to bottom. The weakest gray level shows the sample mask. 



close to the standard 'concordance' power spectrum), gen- 
erated the realization on a grid with the step of 1 Mpc//i, 
and smoothed the field with a Gaussian of a = 2 Mpc/h. 
The original density distribution, the scaling distributions 
and the wavelets are shown in Figs. [7| and [H] for a slice at 
z = 34 Mpc/h, at about the middle of the sample volume. 

The Minkowski functional Vz for the wavelets is shown 
in Fig. [5] In order not to overcrowd the figure, we do not 
show the theoretical predictions. The functional has been 
rescaled to show all functionals together in a single diagram 
and scaling factors are shown in the labels. 

As we see, the lower the wavelet order, the larger values 
of V3 are obtained (this is also true for the other function- 
als); this is expected, as higher orders represent increasingly 
smoother details of the field. The values of the functionals 
for the zero-order wavelet are always higher than those for 
the full field, as it includes only the high-resolution details 
that the iso-levels have to follow. Also, we have seen that 
the lower the order of the functional, the smaller are the 
distortions from Gaussianity. 

The distortions of the third Minkowski functional are 
the largest. The functional for the zeroth-order wavelet 
(curve WO in Fig.|5J shows argument compression, the result 



of insufficient spatial resolution of the smallest details of the 
field. The functional for the first-order wavelet is close to 
Gaussian, as is the functional for the total realization, but 
the second-order curve 8W2 shows strong distortions, due 
to a small number of independent resolution elements in the 
volume, and to the small height of the slice. The character- 
istic volume of these elements is (4 x 2 2 ) 3 = 4096 Mpc 3 //i 3 , 
and their number is about 670. The functional for the third 
wavelet order is already completely dominated by noise. 

So, we can estimate Minkowski functionals with a high 
precision, and the border corrections work well. The most 
difficult part at the moment is the scale separation in ob- 
served samples. The sample geometries are yet slice-like, 
limiting the range of useful scales by the mean thickness 
of the slice. Complex sample borders are also a nuisance 
when applying the wavelet cascade. In order to take ac- 
count of these difficulties, we have yet to resort to running 
Monte-Carlo simulations of Gaussian realizations of a right 
power spectrum, and to compare the obtained distributions 
of the functionals with the functionals for the galaxy data. 
We hope that for the future surveys (e.g., the full SDSS), the 
data volume will be large enough to do without Monte-Carlo 
runs. 
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Figure 11. Multi-scale decomposition of the 2dfl9S volume-limited sample, for the z = 97 Mpc/h slice, and for the grid step v2 Mpc/h. 
The original density is shown at top right, the last scaling order at bottom left, the wavelet orders in between. The wavelet orders increase 
from right to left and from top to bottom. The weakest gray level shows the sample mask (in most cases; sometimes the mask is missing 
and sometimes the gray level is wider than the mask). 



5 MORPHOLOGY OF THE 2DF19 SAMPLE 

Having developed all necessary tools, we apply them to the 
2dF19 volume-limited sample, separately for the NGC and 
SGC regions. 

We show selected slices for the two subsamples, first 
(Figs. EHUD - The Northern slice was chosen to show the 
richest super-clust er in the 2dfGRS N GC, super-cluster 126 
(middle and low, l|Einast,o et, alJlTo^) 1 ). The Southern slice 
has the maximum area in the z = const slices of this sample. 
We choose the gray levels to show also the mask area. 

Figs. 11211151 show the results - the Minkowski func- 
tionals for the wavelet decompositions of the 2dFN19 and 
2dFS19 subsamples. 

First, we show the summary figures: the density V3 for 
the original data and for wavelet orders from to 3, and for 
grid unit of y2 Mpc/h. Higher wavelet orders are not usable 
- there the boundary rules used for the wavelet cascade in- 
fluence strongly the results, and the number of independent 
resolution elements becomes very small, letting the func- 
tionals to be dominated by noise. The wavelet orders used 
span the scale range 2-22.6 Mpc/h. In order to show all the 
densities in the same plot, we use the mapping 



logn(v 3 ) = sgn(« 3 )log(l + |a>3|)- 

It is almost linear for | V3 1 < 1, logarithmic for \v-j\ > 1, 
and can be applied to negative arguments, also. The density 
in Fig. I12l is shown with dotted lines. For reference, the two 
thick lines show the Gaussian predictions for small and large 
functional density amplitudes for the logn(w3) mapping. The 
first glance at the figures of the fourth Minkowski functional 
reveals that none of the wavelet scales shows Gaussian be- 
havior. 

In order to estimate the spread in the values of the 
functionals, we ran about 100 Gaussian realizations for ev- 
ery sample and grid, generated wavelets and found the 
Minkowski functionals. The power spectrum for these re- 
alizat ions was chos e n as described in Sec. 0] above (see 
iKlypin fc Holtzmanl jl997T) ). and smoothed by a Gaussian 
of a = 1 (in grid units). This is practically equivalent to the 
B3 extirpolation used to generate the observed density on 
the grid. Now, if the observational MF-s lie outside the limit- 
ing values of these realizations, we can say that the Gaussian 
hypothesis is rejected with the p- value less than 1%. We also 
calculated the multiscale functionals for a s et of 22 mock 
samples, specially created for the 2dFGRS iNorberg et alJ 
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Figure 12. Summary of the densities of the fourth MF 1)3 for the data and all wavelet orders for the 2dFN19 sample (left) and 2dFS19 
sample (right), in the logn mapping — the higher the wavelet order (indicated by labels), the lower the density amplitude. Thick lines 
show reference Gaussian predictions. Thin lines stand for the a/2 Mpc//i grid. 



(2002). The mock catalogs were extracted from the Virgo 
Consortium ACDM Hu bble volume simu lation, and a bias- 
ing scheme described in ICole et alj iTflSa was used to pop- 
ulate the dark matter distribution with galaxies. 

Fig. 1131 and 1141 show respectively the densities of the 
fourth Minkowski functional for the Northern and the South- 
ern Galactic caps. The MF density v$ corresponding to the 
data is plotted with a continuous line while error bars cor- 
respond to the total variation for mocks. The minimum and 
maximum limits for Gaussian realizations are plotted with 
dotted lines. 

The Gaussian realizations show very small spread, and 
are clearly different from the U3 Minkowski Functional of the 
observational samples. Results from mocks are much closer 
to the data. 

Fig. H3l shows clear non-Gaussianity for the north galac- 
tic cap at a high confidence level. Gaussian realizations are 
not much deformed by the combination of boundary condi- 
tions (wavelets) and border corrections (functionals) effects. 
In this figure we can appreciate that with respect to this 
MF, mocks follow data well for smaller density iso-levels, 
but deviate around v = 1 ; we have seen similar effects be- 
fore l|Martmez et alJboOfl ). An interesting detail is the knee 
around v — 0.5, seen both in the data and in the mocks, but 
not for Gaussian realizations. The latter fact tells us that it 
is not caused by the specific geometry of the data sample. 
The grid step was lMpc/ft. 

Fig. 1141 for the Southern data shows also clear non- 
Gaussianity, but no strong features like the northern one 
(it describes larger scales, as this wavelet is based on the y/2 
Mpc/h step grid). An interesting point is that mocks follow 
the data curve almost perfectly here, much better than for 
the Northern sample. 

Fig. 1151 shows the V3 MF density at the first wavelet 
scale for both north and south slices. Again, it is clearly non- 
Gaussian. However, it is remarkable how well the functionals 
for both volumes coincide. This shows that the border and 
boundary effects are small, and we are seeing real features 
in the density distribution that we have to explain. Also, 
the mocks follow the data rather well. This means that the 
structure (and galaxy) formation recipes used to build the 
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Figure 13. The density of the fourth MF V3 for the wavelet order 
2 for the 2dFN19 sample (full line). Dotted lines show the minima 
and maxima of 102 Gaussian realizations, and bars show the full 
variation in a sample of 22 mock catalogues. 
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Figure 14. The density of the fourth MF V3 for the wavelet order 
2 (grid y/2) for the 2dFS19 sample (full line). Dotted lines show 
the minima and maxima of 108 Gaussian realizations, and bars 
show the full variation in a sample of 22 mock catalogues. 

mocks already implicitly include mechanisms responsible for 
these features. 

It is useful to compare our results with the recent care- 
ful analysis of the topology of the SDSS galaxy distribu- 
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Figure 15. The densities of the fourth MF V3 for the wavelet 
order 1 for the 2dFS19 and 2dFN19 samples (full lines). Bars 
show the full variation in two samples of 22 mock catalogues. 

tion by |Park"et al. (2005f). They used Gaussian kernels to 
find the density distribution, and as a result, their genus 
curves (see, e.g., their Figs. 6 and 8) are close to those 
of Gaussian fields. They describe deviations from Gaus- 
sianity by moments of the genus curve, taken in carefully 
chosen v intervals, and normalized by corresponding Gaus- 
sian values. Our Minkowski functionals differ so much from 
Gaussian templates (see Fig. I12H that we cannot fit a refer- 
ence Gaussian curve. The only analogue we can find is the 
shift of the genus curve Av, that we estimate by fitting an 
expression V3(v) = A {{u — Au) 2 — l) exp \ — {y — Av) 2 /2) 
to our results. The values of the shift corroborate the vi- 
sual impressi on of strong diffe rences between our results 
and those of iPark et al.l (120051) . While they find that Av 
lies in the interval [—0.1,0.26] (for the scale range Rg £ 
[4.5, 11.0] Mpc/ft, their table 2), our Av assumes values be- 
tween — 1.3 and —0.2, for approximately the same scale in- 
terval (A € [4, 22.6] Mpc/fe). As the morphology of the 2dF- 
GRS and the SDSS should not differ much, the difference 
is clearly caused by different kernels, Gaussians compared 
to compact wavelets. Dropping the conventional Gaussian 
kernels makes the discriminative force of the morphological 
tests considerably stronger. 



6 CONCLUSIONS 

The main results of this paper are: 

(i) We have shown how to compute the MF, taking into 
account both the biases due to the discrete grid from the 
Crofton method and the border effects related to complex 
observational sample volumes. 

(ii) Our experiments have shown that the multiscale MF 
functionals of a Gaussian Random Field have always a Gaus- 
sian behavior, even in case where the field lies within com- 
plex boundaries. Therefore, we have established a solid base 
for calculating the Minkowski functionals for real data sets 
and their multi-scale decompositions. 

(iii) We found that both the observed galaxy density 
fields and mocks show clear non-Gaussian features of the 
morphological descriptors over the whole scale range we have 
considered. For smaller scales, this non-Gaussianity of the 
present cosmological fields should be expected, but it has 
been an elusive quality, not detected in most of previous 



papers (see, e.g., iHovle et al.l (I2002I) . IPark et al.l J2005M . 

But even for the largest scales that the data allows us to 
study (about 20 Mpc/h), the density fields are yet not Gaus- 
sian. We believe that the Gaussianity reported in the papers 
cited above could be just a consequence of oversmoothing 
the da ta. This effect was clearly described in lMartmez et all 
J2005I) . 

(iv) The mocks that are generated from initial Gaus- 
sian density perturbations by gravitational evolution and 
by applying semi-analytic galaxy formation recipes, are 
pretty close to the da ta. However, as in a previous study 
llMartfnez et alJl2005l) . we confirm a discrepancy around 
v = 1 between the mocks and the data. This analysis clearly 
shows that there are more faint structures in the data than 
in the mocks, and clusters in the mocks have a larger inten- 
sity than in the real data. 
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APPENDIX A: THE A TROUS ALGORITHM 
AND GALAXY CATALOGUES 

A good description of the a trous algorithm and of its ap- 
plications to image proce ssing in astronomy can be found in 
IStarck & Murt agh (2002). Readers intereste d in the math - 
ema tical ba s is of the algorithm can consult iMallad (ll99St) 
and IShensal (Il992f) . We give below a short summary of the 
algorithm and describe the additional intricacies that arise 
when the algorithm is applied to galaxy catalogues (point 
data) . 

We start with forming the initial density distribution 
do on a grid. In order to form the discrete distribution, we 
have to weight the poi nt data (extirp olate), using the scaling 
kernel for the wavelet (fMallat 1999). As we shall use the B 3 
box spline as the scaling kernel, the extirpolation step is: 



d (0 V) 



/9(x)_Bg (x — rij)d x, 



(11) 



where rii = (n)i — (xi,yi,Zi) is a grid vertex, p(x) is 
the original density, delta-valued at galaxy positions, and 
-Bg^x) is the direct product of three S3 splines: 



B< 3) (x) = B 3 (x)B 3 (y)B 3 (z), 

where (x) — (x,y,z). The B 3 spline is given by 



■41a;- 



1| 3 +6|a;| 3 



-4|a; + l| 3 + |: C + 2| 3 l 



As this function is zero outside the cube [—2, 2] 3 , every data 
point contributes only to its immediate grid neighbourhood, 
and extirpolation is fast. 

The main computation cycle starts now by convoluting 
the data d with a specially chosen discrete filter h(u\ : 



d' 



Eh 
a ( fc ) a (n) + 2'(fc)- 

(*0 



(12) 



Here I stands for the convolution order (octave), the 3- 
dimensional filter = hih m h n , (k) = (l,m,n) is the di- 
rect product of three one-dimensional filters hi — {1/16, 1/4, 
3/8, 1/4, 1/16}, for i g [-2,2]. Two points should be noted: 

• As the filter is the direct product of the one-dimensional 
filters, the convolution can be applied consecutively for each 
coordinate, and can be done in place, with extra memory 
only for a data line. 

• The data index (n)+2 I (k) in the convolution shows that 
the data is assessed from consecutively larger regions for 
further octaves, leaving intermediate grid vertices unused. 
This is equivalent to inserting zeroes in the filter for these 
points, and this is where the name of the method comes 
from (a trous is "with holes" in French). This makes the 
convolution very fast, as the number of operations does not 
increase when the filter width increases. 

The filter hi is satisfies the dilation equation 

i*(f)-£- 



h k B 3 (x — k). 



After we have performed the convolution 112L we find 
the wavelet coefficients ti/ J ' for the octave J by simple sub- 
straction: 



.,( J ) 



d 



(J) 



d 



(J+i) 



"(n) - "(n) - "(n) ■ (I 3 ) 

The combination of steps (1121131 is equivalent to convolution 
of the data with the associated wavelet ip^\x), where 

tl>(x) = 2B 3 (2x)- B 3 (x). (14) 

Repeating the sequence 112I13H we find wavelet coefficients 
for a sequence of octaves. The number of octaves is, evi- 
dently, limited by the grid size, and in real applications by 
the geometry of the sample. 

We have illustrated the wavelet cascade in the main 
text by application to real galaxy samples. As our wavelet 
amplitudes were obtained by subtraction, we can easily re- 
construct the initial density: 



.CO 
'(«)• 



(15) 



(n) (n) 

3=0 

Here the upper indices show the octave, the lower indices 
denote grid vertices, and d^ J+1 ^ is the result of the last con- 
volution. This formula can also be interpreted as the decom- 
position of the original data (density field) into contributions 
from different scales - the wavelet octaves describe contri- 
butions from a limited (dyadic) range of scales. 

Here we have to note that while the scaling kernel 
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Figure 16. The square of the Fourier transform of the wavelet 
for two neighbouring octaves. 



$(x,y,z) = B 3 (x)B 3 (y)B 3 (z) 

is a direct product of three one-dimensional functions, it 
is surprisingly almost isotropic. Its innermost iso-levels are 
slightly concave, and outer iso-levels tend to be cubic, but 
this happens at very low function values. In order to char- 
acterise the deviation from anisotropy, let us first define the 
angle- averaged scaling kernel 

*(r) = f {r)${r,0,4>)dS, 
J s 

where S(r) is a spherical surface of a radius r. The 
anisotropy can now be calculated as the integral of the abso- 
lute value of the difference between the kernel and its angle- 
averaged value: 



$(x,y,z) — $(\/x 2 + y 2 + z 2 ) 



dxdydz = 0.030. 



As the integral of the kernel itself is unity, the deviation is 
only a couple of per cent. 

A similar integral over the wavelet profile gives the value 
0.052. Here the natural scale, the integral of the square of the 
wavelet profile, is also unity, so the deviation from isotropy 
is small. 

Isotropy of the wavelet is essential, if we want to be sure 
that our results do not depend on the orientation of the grid. 
This is usually assumed, but with a different choice of the 
scaling kernel this could easily happen. 

As our wavelet transform is not orthogonal, there re- 
main correlations between wavelet amplitudes of different 
octaves. The Fourier transform of the B 3 scaling function is 



B 3 (u) = 



I sin(u/2) 



V w / 2 

The Fourier transform of the associated wavelet (1131 is 

w(lu) = B 3 (u/2) -B 3 (u). 

We show the square of the Fourier transform of the 
wavelet for two neighbouring octaves in Fig. 1161 As we see, 
the overlap between the octaves is not large, but substantial. 
This is the price we pay for keeping the wavelet transform 
shift invariant. We can now compare wavelet amplitudes for 
different octaves (scale ranges) at any grid vertex, but we 
have to keep in mind that the separation of scales is not 
complete. It may seem an unpleasant restriction that the 




Figure 17. Relative errors of border-corrected densities of the 
second MF v\ for a realization of a Gaussian random field in the 
2dfl9N sample mask. The case of the border correction chain is 
shown by solid line, the 'raw' correction case - by dotted line. 



wavelet scales have to increase in dyadic steps. It is not, in 
fact, as one can choose the starting scale (the step of the 
grid) at will. 

Before applying the wavelet transform to the data, we 
have to decide how to calculate the convolution 1121 near the 
spatial boundaries of the sample. Exact convolution can be 
carried out only for periodic test data, and spatially limited 
data need special consideration. For density estimation, a 
useful method to deal with boundaries is to renormalise the 
kernel. This cannot be done here, as renormalisation would 
destroy the wavelet nature of our convolution cascade. The 
only assumptions that can be used are those about the be- 
haviour of the density outside the boundaries of the sample. 
Let us consider, for example, the one-dimensional case and 
the data d(i) known only for the grid indices i > 0. The 
possible boundary conditions are, then: 



f ° 

d(i;i<0) = < d(0) 

1 d(-i) 



zero boundary 
constant boundary 
reflecting boundary. 



The constant boundary condition is rarely used; the most 
popular case seems to be the reflecting boundary. For brick- 
type sample geometries, where the coordinate lines are per- 
pendicular to the sample boundary, this condition gives good 
results. However, in our case the sample boundary has a 
complex geometry, and reflections from nearby boundary 
surface details would soon interfere with each other. 



APPENDIX B. COMPARING BORDER 
CORRECTIONS 

We noted above (Sec. 12.41 that alongside with the border 
correction chain 1101 there is another possibility to correct 
for borders. In this case we ignore the vertices in the mask 
and do not build any basic elements if one of the vertices 
belong to the mask. This also means that we do not have to 
build the basic elements in the mask region. The latter fact 
makes the algorithm faster (about twice faster for the 2dF 
data), as the data region occupies usually only a fraction of 
the encompassing brick). 

We call this method the 'raw' border correction and 
compare it with the border-correction algorithm 1101 on the 
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Figure 18. Relative errors of border-corrected densities of the 
third MF V2 for a realization of a Gaussian random field in the 
2dfl9N sample mask. The case of the border correction chain is 
shown by solid line, the 'raw' correction case - by dotted line. 



than that for the 'raw' border correction. Thus we use this 
chain throughout the paper. 

It is useful to recall, though, that the border correc- 
tions (llOt are based on the assumption of homogeneity and 
isotropy of the data, which may not always be the case. The 
'raw' border corrections do not rely on any assumptions, and 
are therefore useful for verifying the results obtained by the 
correction chain. 
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Figure 19. Relative errors of border-corrected densities of the 
fourth MF v-j for a realization of a Gaussian random field in the 
2dfl9N sample mask. The case of the border correction chain is 
shown by solid line, the 'raw' correction case - by dotted line. 



example of the realization of a Gaussian random field for the 
2dF NGC region, smoothed by a Gaussian of a = 3 Mpc//i 
to ensure that we resolve the density distribution. (We used 
the same realization to compare the border-corrected and 
uncorrected case, in Sec. 12.41 1 We combine the encompass- 
ing Gaussian brick with the 2dFN19 mask, calculate the 
Minkowski functionals for both border correction methods, 
and compare them with the functionals found for the peri- 
odic brick. We show below the relative errors of the func- 
tionals, defined as 



max„ |v? 



where (u) are the densities of the functionals for the brick. 
We cannot use v\(v) themselves to normalise the errors, as 
their values pass through zero, so we use their maximum 
absolute values. 

The relative differences are shown in Figs. 1171191 We 
see that both correcting methods give the results that do 
not differ from the true densities more than 10% (12% for 
V3). We see also that the border correction chain [1 1 ( H gives 
always better estimates of the functionals; the maximum 
error is 3-4%, and the error is about three times smaller 



